Importance of including small body spin effects in the modelling of extreme and 

intermediate mass-ratio inspirals 



E.A. Huerta 1,2 and Jonathan R. Gair 1 

1 Institute of Astronomy, Madingley Road, CBS OH A Cambridge, UK 
2 Department of Physics, Syracuse University, Syracuse, NY 13244> USA. 
(Dated: January 21, 2013) 

In this paper we explore the ability of future low-frequency gravitational wave detectors to mea- 
sure the spin of stellar mass and intermediate mass black holes (IMBHs) that inspiral into spinning 
supermassive black holes (SMBHs). We describe a kludge waveform model based on the equa- 
tions of motion derived by Saijo et al. [l[ for spinning BH binaries, augmented with spin-orbit 
and spin-spin couplings taken from perturbative and post-Newtonian (PN) calculations, and the 
associated conservative self-force corrections, derived by comparison to PN results. We model the 
inspiral phase using accurate fluxes which include perturbative corrections for the spin of the in- 
spiralling body, spin-spin couplings and higher-order fits to solutions of the Teukolsky equation. 
We present results of Monte Carlo simulations of parameter estimation errors, computed using the 
Fisher Matrix formalism, and also of the model errors that arise when we omit conservative cor- 
rections from the waveform template. We present results for the inspirals of spinning BHs with 
masses = 10M Q , 1O 2 M , 10 3 M Q , 5 x 10 3 Af e , into SMBHs of specific spin parameter q = 0.9, 
and mass M = 10 6 Mq. The analysis shows that for intermediate-mass-ratio inspirals, e.g., a source 
5 x 10 3 M© + 1O 6 M0 observed with a signal-to-noise ratio (SNR) of 1000, observations made with 
LISA will be able to determine the inspiralling BH mass, the central SMBH mass, the SMBH spin 
magnitude, and the magnitude of the spin of the inspiralling BH to within fractional errors of 
~ 10 -3 , 10 -3 , 10 -4 , 10%, respectively. We also expect to determine the location of the source in 
the sky and the SMBH spin orientation to within ~ 10~ 4 steradians. However, for extreme-mass- 
ratio inspirals, e.g., a WMq + 10 6 Mq system observed with SNR of 30, LISA will not be able to 
determine the spin magnitude of the inspiralling object, although the measurement of the other 
waveform parameters is not significantly degraded by the presence of spin. We also show that the 
model errors which arise from ignoring conservative corrections become significant for mass-ratios 
above ~ 10~ 4 , but including these corrections up to second PN order may be sufficient to reduce 
these systematic errors to an acceptable level. 
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I. INTRODUCTION 

Future low-frequency gravitational wave (GW) detectors, such as LISA [2|, will provide information about the 
astrophysical properties of a great variety of sources, including merging SMBHs, stellar mass compact object (CO) 
binaries in our galaxy and the inspirals of stellar mass COs into SMBHs in the centres of galaxies. The scientific payoffs 
of such observations will shed light on the spin and mass distributions of SMBHs, the mass distributions of inspiralling 
stellar-mass objects in galactic centres and the cosmological rate of mergers, and unprecedented information on the 
properties and evolution of BHs in the Universe. 

In addition to providing this unique new information about astrophysical objects in the Universe, detectors like 
LISA will enable us to carry out tests of general relativity in the strong-field regime close to the horizon of massive 
BHs. These tests will allow us to determine whether the geometry of the BH is the Kerr metric of relativity, as we 
expect, or if it is some other kinds of massive body, such as a boson star, naked singularity, soliton star, etc. The 
extreme-mass-ratio inspirals (EMRIs) of stellar-mass compact objects into much more massive SMBHs in the centres 
of galaxies are ideally suited for such tests, as many hundreds of thousands of waveform cycles are generated while 
the smaller body is orbiting in the strong field regime close to the SMBH. 

Ryan Q showed, for nearly-circular nearly-equatorial EMRIs, that the phase of the GWs observed will provide a map 
of the spacetime structure outside the central BH. Such a map might also be extracted for eccentric-inclined EMRIs 
[H, but with a few additional complications, since the eccentricity and inclination cannot be directly determined, but 
only inferred from the observed frequencies Q. These maps can be compared to the Kerr metric in order to detect 
deviations from the general relativistic predictions. If deviations are observed, these could arise from a difference in 
the nature of the central object or a failure in the theory of relativity, but they might also have an astrophysical 
origin, e.g., the presence of gas in the system. However, most massive holes are in galaxies with non-active nuclei, 
and will be surrounded by only a tenuous disk with advection-dominated accretion flow. Narayan Q has shown 
that in such environments accretion drag should be totally negligible, and hence the orbital evolution will be entirely 
driven by radiation reaction. Massive BHs in active galactic nuclei will be surrounded by thin accretion disks. In that 
case, accretion-disk drag will be significant, which could leave an imprint on the EMRI signal, but in general such 
astrophysical effects are qualitatively different from the effect of a deviation in the spacetime structure, e.g., they can 
drive increasing eccentricity Q- 

If we are to detect small differences in the gravitational wave phasing, it is necessary to develop accurate waveform 
templates for the pure Kerr-inspiral model, against which to test the observations. Great progress has been made in 
the modelling of EMRI systems and the understanding gained from this work has driven the development of algorithms 
for EMRI detection and data analysis [H-[l3j]. The formalism to compute accurate signals is already well developed. 
Teukolsky [3] and Sasaki et al. [l5[ developed a mathematical framework to study first-order perturbations of Kerr 
BHs. However, this approach is computationally expensive and has left a need for the development of computationally 
inexpensive waveform templates that capture the main features of the true signals, and that can be used for source 
detection in data analysis algorithms. In this paper we will use one particular approach, the so-called "numerical 
kludge" waveform model. This approach combines a prescription to capture the radiative dynamics of Teukolsky- 
based evolutions (l6| . using accurate expressions for the fluxes of energy and angular momentum to compute the 
inspiral trajectories with an approximate flat-spacetime model to construct the waveform from the trajectory [l7j . It 
has been shown that this approach gives waveforms whose overlap with fully accurate general relativistic waveforms 
is greater than 95% over a large portion of the parameter space [16]. This approach has recently been improved 
for the case of circular-equatorial EMRIs, by including conservative self-force corrections at 2PN order [l8|]. These 
corrections are important because they affect the phasing of the waveform over time, and neglecting them may lead 
to several cycles of discrepancy in the waveform template over the course of an inspiral. The self-force program has 
not yet obtained general relativistic conservative self-force corrections for Kerr inspirals, so in [l8[ these corrections 
were derived by comparison to PN results and requiring that the "numerical kludge" model reproduced PN results in 
the weak field. Using this improved waveform model, we showed that the inclusion of such corrections in the model 
was not essential for source detection, but we may have to include them for accurate parameter determination for 
EMRIs. In this paper, we consider another effect, namely the influence of the spin of the inspiralling body on Kerr 
circular-equatorial inspirals and assess the importance of including this parameter in EMRI models. 

This analysis is important because it is expected that most astrophysical BHs will have significant spin and therefore 
the true gravitational waveforms will include small body spin effects. It has been argued that such effects may not 

ta significant role for signal detection and parameter estimation for EMRIs with mass ratios r\ = /J./M ~ 10 -5 
This is because the spin of the inspiralling object is suppressed by a factor of r\ in the equations of motion 
and hence the inclusion of small body spin effects will modify the orbital evolution by at most a few radians over a 
year. Studies on the effect of the self-force on spinning objects that inspiral into Schwarzschild BHs have shown that 
the accurate determination of the local self-force may not be needed for a determination of the orbital evolution of 
EMRIs, and its omission will only introduce a small error in the determination of the spin rate of the companion [20| . 
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However, the spin effects become more important as the mass ratio, 77, is increased and can therefore be significant 
for intermediate mass ratio inspirals (IMRIs) with 77 ~ 10~ 3 (e.g., the inspiral of a ~ 10 3 M o intermediate mass BH 
(IMBH) into a ~ 10 6 M Q SMBH). Even if small body spin effects are only marginally important for detection, their 
inclusion in waveform templates may be important for parameter estimation. We also need to determine at what 
mass ratio small body spin effects become important for signal detection, i.e., when GW observations will be able to 
measure the spin of the inspiralling body. We address this question here by developing a "numerical kludge" waveform 
model that includes small body spin effects using the mathematical machinery developed by Saijo et al. and by 
including conservative self-force corrections for spin-orbit and spin-spin couplings using the same method that was 
employed to perturbative conservative corrections in (l8| . We compute, for the first time, the mass-ratio threshold 77 
above which low-frequency gravitational wave detectors like LISA will be able to accurately measure the spin of a BH 
inspiralling into a SMBH. We will show that future observations will not only provide an unprecedented census of the 
spin distribution of SMBHs, but could also yield an accurate census of the spin and mass distributions of intermediate 
mass BHs with mass fj, ;> 10 3 M Q . 

This paper is organized as follows. In Section [H] we introduce the waveform model used to study the inspiral of 
spinning BHs into SMBHs. We construct kludge asymptotic observables to implement the radiative and conservative 
components of the self-force, including spin-orbit and spin-spin couplings, by comparison to PN results. Section IlIII 
presents the wave generation prescription used in our studies, along with the implementation of LISA's response 
function and Doppler phase modulations in our waveform model. In Section ITVl we review the basic elements of signal 
analysis, which we use in Section [V] to assess the accuracy with which LISA observations might be able to determine 
the parameters of a representative sample of binary systems. In Section [VII we describe the formalism developed by 
Cutler & Vallisneri [2l| to estimate the theoretical or "model" errors that would arise from omitting conservative 
corrections in the waveform template. We summarize our conclusions in Section IVIII 



II. KLUDGE WAVEFORM WITH SMALL BODY SPIN EFFECTS 



The numerical kludge waveform model developed by Babak et al. |16| has been very successful as a model for 
the GW emission from EMRIs, as it is able to capture the main features of the inspiral waveform in the strong field 
regime. In this scheme, one assumes that the inspiralling object instantaneously follows a Kerr space-time geodesic, 
and slowly evolve the parameters of the geodesic using a prescription for the rate of change of energy, angular 
momentum and Carter constant that is based on PN expressions, augmented by fits to Teukolsky-based evolutions 
(l7| . This results in a trajectory for the inspiralling body in Boyer-Lindquist coordinates. The gravitational waveform 
is then constructed by identifying these coordinates with spherical polar coordinates in a pseudo-flat space and using 
a weak-field emission formula. This model has already been used to explore the accuracy with which LISA will be able 
to measure the parameters of Kerr circular equatorial EMRIs [l8|, an d to estimate the importance of conservative 
and radiative self-force corrections on parameter estimation and detection. In the kludge models to date, the spin of 
the stellar mass CO was ignored. For the purpose of providing a better theoretical template, and to find out under 
which circumstances the spin of the CO can be measured, in this paper we will develop a new kludge model that 
incorporates this additional parameter. 

To achieve this, we will augment the standard "numerical kludge" model using the equations of motion of a spinning 
particle in the equatorial plane (0 = n/2) of a Kerr BH, as derived by Saijo et al. We will include small body spin 
corrections of two different natures: i) first-order conservative corrections to amend the orbital phase evolution, and; 
ii) second order radiative corrections in the fluxes of energy and angular momentum to evolve the geodesic parameters 
of the inspiralling object. 

For a particle with spin angular momentum Si = s/iz, aligned with the central Kerr BH spin (S2 = aMz) and the 
orbital angular momentum L z , the spin vectors remain constant and the equations of motion take a form similar to 
the Kerr geodesic equations, namely, 
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and (t 7 p,9,ip) are Boyer-Lindquist coordinates, r is the particle's proper time, and E = E/^j, and Jj = J z /^ are the 
conserved energy and total angular momentum, respectively (see Eqs. (2.10) of [![). The radial motion of the spinning 
particle can be understood by re-writing the function R s in the form 



R s = B(p)[E - Exip, J Z )][E - E 2 {p, J z )], 
where the roots E± t 2 of R s = are found by solving 
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The effective potential is defined as the minimum allowed value of the particle energy at radius p and so Eq. 
implies, 



(8) 



Note that we have taken the positive square root because the motion of the spinning particle should be future directed. 

We can now derive the expressions for the energy and angular momentum of a spinning particle, following a circular 
equatorial orbit by ensuring that the function R s (p), and its radial gradient, R s (p), vanish at the same point. 

The equations of motion derived by Saijo et al. [lj are valid only to linear order in the spin of the small body, so we 
shall quote the expressions for the energy and angular momentum of the spinning particle at the same order, namely, 



E _ r 2 - 2r ± (q + s/r)^r + iqs/r - 5qs/2r 
^ r^J r 2 — 3r ± (2q + 3s/r)y?r + Sqijr — Qqijr 
L z _± v / r + 3qs/r (r 2 + q 2 + qs(r + l)/r) - 2rq + sr(r - |) 
r^J r 2 — 3r ± (2q + 3s/r)y^r + Sqijr — 6qs~Jr 

where r — p/M , s — s/M — r]x, with x the dimensionless spin parameter of the inspiralling BH, r\ = fi/M is the mass 
ratio as before, and q = a/M, with a the spin of the central BH. 

Furthermore, for Schwarzschild BHs, q = 0, we find that at linear order in s, 
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Eq. dXDJ corrects a typo in Eq. (B18) of [22J. 

To obtain the orbital evolution of the CO, we need to calculate the evolution of the energy E and angular momentum 
L z . We can evaluate these quantities by equating their rate of change with the flux carried away by the GWs, E and 
L z . We will use the radiation fluxes derived by Gair & Glampedakis [ItJ augmented with accurate BH perturbation 
theory (BHPT) results that include small body spin corrections [23j], i.e., 



• 32^ 2 /l\ 5 f 1247 AN / 73 25 \ /l\ 3/2 ( 44711 33 , 71 \/l x2 

+ higher order Teukolsky fits J , 

■ 32u 2 {l\ r/2 ( 1247 /l\ / 61 19 \ /1\ 3/2 / 44711 33 , 59 \/l^ 2 
L ^-^M{r) ^-WUj + r-^-T^JUj + (-^072- + 16 9 + Y™ 

+ higher order Teukolsky fits j> . (12) 

The "higher order Teukolsky fits" are given in [l7|. We do not give these explicitly here, as they are not needed to 
derive the conservative corrections. However, we will include them to evolve orbits when we generate waveforms. 

To evolve a circular orbit for a spinless particle in a Kerr background, we need only to specify the angular momentum 
or the energy flux, as they are related by the 'circular goes to circular' rule [13] 

E(r, X -> 0) = ± 1 L z (r, x -> 0) = Q(r, X -> 0)L,(r, X -> 0), (13) 
r J / z ± g 

where d0/d£ = r2(r), is the azimuthal velocity of the orbit. We can use exactly the same scheme for spinning particles, 
as Tanaka et al. [231 ] showed that the assumption that a circular orbit remains circular under radiation reaction is 
consistent with the energy and angular momentum loss rates at linear order in the spin of the particle. Thus, Eq. 
(fT3| now reads [23J, l25( 

*< r »- ± ?^i( 1 -5«?^) t -w- <14) 

Furthermore, the evolution in time of the radial coordinate is given by 

dr • dr • 

'-w'-ac*- (15) 

In the non-spinning case, \ = 0, using the exact geodesic expression for dr/dL z generates inspirals that are closer to 
Teukolsky based evolutions than expanding the above expression at 2PN order [181 ] and we will adopt this approach 
here as well. However, we will need the 2PN expression in the following, which is 

dr 64 r) ( d\ 3 f 743 fl\ ( 133 35?7x\/l x3/2 
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Up to this point our analysis has been incomplete as we have only considered the three shape integrals of the motion. 
There are also three positional constants of the motion, which label the position of the test particle along the geodesic 
trajectory at some fiducial time. The self-interaction between the inspiralling body and the central BH has two main 
pieces. So far we have considered only the dissipative or radiative self-force, which drives the evolution of the shape 
constants. However, there is also a second piece, the conservative self-force, which affects the positional constants 
of the motion. Effectively, this changes the frequency of an orbit at a given radius but does not cause the orbit to 
evolve. 

The conservative self-force has two parts: one is oscillatory and averages to zero, but there is also a secular piece 
which leads to accumulation of a phase error over time. In our model, we will include this effect by amending the 
evolution equation for the <f> frequency as follows, 



'd(f> 
dt \~dt 



gco 



1 + SQ . (17) 



This equation includes the phase derivative for a geodesic, labeled by the subscript "geo" which is given implicitly 
in Eq. (|14l) . and a frequency shift which will depend on the instantaneous orbital parameters. A problem arises 
here because to compute the necessary frequency shifts within our framework, i.e., BHPT, would require self-force 
calculations. At present, Barack & Sago have computed gravitational self-force corrections for particles moving on 
circular and eccentric bound geodesic orbits around a Schwarzschild BH in the Lorenz gauge (2(| [27| • Warburton & 
Barack have recently computed the self-force on a scalar charge for Kerr circular and eccentric equatorial orbits in 
the same gauge [28l . [29j . The extension of this work to Kerr eccentric inclined geodesic orbits is still a challenging 
endeavour mainly because there is not a formal framework to deal with Lorenz-gauge metric perturbations in the 
frequency domain — the natural arena in which these calculations are carried out. 

However, we do know conservative corrections in the PN framework up to 2PN order which include spin-orbit, 
spin-spin couplings and finite mass contributions [3Qt | . We can combine these expressions with the radiative self-force 
obtained by Tanaka et al. [2j| based on the Teukolsky and Sasaki-Nakamura formalisms for perturbations around 
a Kerr BH, which includes terms of order q 2 , r]x, qvX- Using these results, we will extend the method originally 
proposed by Babak et al. [la ], who computed the 1PN conservative correction for circular orbits in the Schwarzschild 
space-time. This idea was subsequently extended by Huerta & Gair [l^|, who included conservative corrections in 
a kludge model at 2PN order for Kerr circular-equatorial orbits. In this paper we shall extend the latter model by 
including small body spin effects at linear order, and their corresponding conservative corrections at 2PN order. 

The idea is to correct the kludge expressed in a particular coordinate system, in order to ensure that asymptotic 
observables are consistent with PN results in the weak field. In particular, we aim to modify the orbital frequency 
and its first time derivative. By modifying these two quantities, we can both identify coordinates between our kludge 
model and the PN formalism, and find the missing conservative pieces. 

From equation (fT4|) . and using s — r/x, we find that at 2PN order the orbital frequency takes the form 

1 / 1 \ ^ / / 3 \ / 1 \ ^ 3 / 1 \ ^ /l^ ^/ 



which, following Eq. \17\ we now augment by including conservative corrections 



dt M 



3/2 




(19) 

We will use this expansion for J7 geo only to derive the conservative corrections. As with dr/dL z , it has been shown 
that more reliable waveforms can be obtained by including the full geodesic frequency where it is known, see (|17l) , and 
this will be the approach used in section HVl It is inconsistent to include some pieces of the evolution at arbitrary PN 
order, while including the conservative corrections only at 2PN. However, the lower order effects that we are including 
with greater accuracy do have a more significant impact on the waveform [l7| . 
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The expansion in r is an expansion in v 2 = 1/r. To derive the conservative corrections, we choose to leave the time 
derivative of the radial coordinate unchanged and given by equation (|15l) / p^|) . This amounts to a choice of gauge in 
which the i] 2 terms of dr/dt that are not proportional to x vanish. Differentiation of (fl~9| then gives dil/dt for the 
kludge, 



dil 96 7i /l\ 11/2 [ , 1 / 743 /5, 743, 
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(20) 



Note we have used s — r)X, with \ = Si/a* 2 - We may now write down a coordinate transformation to relate our 
coordinates with those used in the PN formalism, namely, 

R l\ fM\ , fM\ 3/2 „ , fM^ 2 
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where R denotes the PN semi-major axis. We can now substitute this expression for the coordinate transformation 
into relations JTSJ and (|2"0"j). 

The final stage of the computation is to compare the expressions for Q and J7, where a dot denotes d/dt, with the 
available PN expansions. The PN expansions are available to higher order in the mass ratio 77, but we keep r\ only to 
the same order as the kludge, Eq. (f2"U|) . The PN expressions for the orbital frequency and its first time derivative are 
given by 30] 



2 \ - 

+ 3r/ \ L-Xi 



where mx = M + m, £ is a unit vector directed along the orbital momentum, x — Xi — Sif/j 2 , q — xi — S2/M 2 . 
Additionally [13], 
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where the constant £ was determined in [j^ . and has the value £ = —243/32. This term guarantees that the PN 
framework and the perturbative approach coincide in the test mass particle limit 77 — > 0. The spin-orbit (3 and 
spin-spin parameters a are given by 
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To be consistent with the perturbative approach outlined above, we shall assume that the spin of the small particle is 
perpendicular to the equatorial plane (this guarantees that the orbit remains circular-equatorial), and parallel to the 
momentum of the central Kerr BH. For prograde orbits, the spin and spin-spin corrections will play a more significant 
role, and so we will focus our analysis on those. We now rewrite these expressions in a convenient way to take the 
small mass-ratio limit, by writing hit = M(l + 77), giving 
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A direct comparison between the expressions for the orbital frequencies and their first time derivatives allow us to solve 
simultaneously for the various coefficients of Eqs. (p~9|) . (|20|) and ([21]) . We find that the non-vanishing parameters are 
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Having found the required corrections, we can now explore whether small body spin effects are important for signal 
detection. We will do this in Section [V] by computing the inspiral trajectory using the relations (|12l) . (IT51) . and using 
the azimuthal frequency amended with conservative corrections, i.e., 



n = 



r 3 / 2 + q 
V ( do + 



2 riX r 2 + q^J I 
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We note that this equation includes corrections up to 2PN order, characterised by the coefficients d etc.. In SectionfVl 
we will assess the importance of these corrections for signal detection and parameter estimation, by turning these 
parameters on and off (i.e., setting a subset of them to zero) in the model. 

Previous studies have ignored small body spin effects, in particular for EMRIs, since they are expected to be small 
for mass rations 77 ~ 10~ 5 . However, even though an accurate measurement of the spin of the small body will be 
unlikely for these systems, we will discuss whether this additional parameter has a bearing on the accuracy with which 
other parameters can be determined. Furthermore, we will study various systems to find the mass ratio threshold 
at which the small spin parameter starts to be measurable through GW observations, and will also explore whether 
small object spin conservative corrections are relevant for parameter estimation and detection. 

The following section briefly outlines the waveform model to be used in these studies. 



9 



III. WAVEFORM MODEL 



Having derived an accurate scheme for the inspiral evolution, we need a prescription for the waveform model. We 
shall generate our waveforms using a flat-space quadrupole wave generation formula applied to the trajectory of the 
inspiralling object in Boyer-Lindquist coordinates, which we identify with spherical-polar coordinates in a flat-space. 
Using this "pseudo-flat space" quadrupole moment approximation is not correct, but this prescription has been seen 
to work well in practice, as shown in (16J. 

In the transverse-traceless (TT) gauge, the metric perturbation hij(t) is given by, 

% (*) = § [PikPjl - \PvPk^j I kl , (29) 

where D is the distance to the source, rjij is the flat Minkowski metric, the projection operator = rjij — hiUj, hi 
is a unit vector in the direction of propagation, and I kl is the second time derivative of the inertia tensor [3l|. In 
the EMRI framework this takes the form I kl = mr % (i)r J (t) , where r l (t) represents the position vector of the CO with 
respect to the MBH in the pseudo-flat space. 



A. Implementation of LISA's response function 



Since most of the SNR of the binary inspirals considered in these studies accumulates at frequencies f <; lOmHz, 
it should be adequate to use the low- frequency approximation to the detector response. Following [32[, the LISA 
response in this regime may be written as 



h a (t) = ^[F+(t)A+(t) + F*(t)A*(t) 



(30) 



where a = I, II refers to the two independent Michelson-like detectors that constitute the LISA response at low 
frequencies. The functions A + ' x (t) are the polarization coefficients given by 



A + = -0+[l + (a • hf], A x = 2a x (a ■ n), 
where a is a unit vector along the SMBH's spin direction, and a + ,a x are given by 

°+4('" 

The antenna pattern functions F^ x are given by 
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F+ = i(l + cos 2 6) sin(20) cos(2V>) + cos6»cos(2<?!)) ^(2^), 



(32) 



F. 



ii 



-(1 + cos 2 6) sin(20) sin(2-0) - cos6»cos(20) cos(2V0- 



(33) 



The various angles in the previous expressions represent the source's sky location in a detector based coordinate 
system, (9,<j>), and the polarization angle of the wavefront, ip. These can be re-written in a fixed, ecliptic-based 
coordinate system. If we denote the source co-latitude and azimuth angles and the direction of a in this fixed 
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coordinate system by (9s, (ps) and (9k,4>k) respectively, then 

1 /3 

cos 9(t) = -cos 9 S - ~y sin S cos[0 o + 2n{t/T) - <j> s ], 

imOs cos[4 + 2w(t/T) - fa] \ 
sm[4> + 2n(t/T) - fa] J' 




2 sin 9s sin[ 
. + 2n(t/T) - cj) K } 



— cos 9(t) [cos 9k cos $s + sin 9k sin 6*s cos(0k — 0s)] J / 

{1 \/3 
- sin 6>k sin 9 S sm.{4> K - <Ps) ^ cos(0 o + 2irt/T) 

{cos 9k sin #s sin <ps — cos 0s sin 9k sin 

\/3 - ] 

— sin(0 o + 2nt/T) (cos 9 S sin 0^ cos 4>k — cos 0^- sin 9s cos 0s) >, (34) 

where 00,^0 are constant angles which represent the orbital and rotational phase of the detector at t — 0. We will 
set both of these to zero in our analysis. Additionally, T is the orbital period, which is 1 year. Barack and Cutler 
[l9j write these expressions in terms of 0l,0l which specify the direction of the CO's orbital angular momentum 
in the ecliptic-based system. In this case the orbits under consideration are circular and equatorial, so the angular 
momentum vector of the orbiting body does not precess about the SMBH spin a, and hence 9k = @l, 4>k — 4>l- 

The last ingredient in the detector response is the Doppler phase modulation. If $(t) denotes the phase of the 
waveform, the inclusion of the Doppler modulation shifts the phase as follows pl| 

$(f) -> $(f) + 2^i?sin 6»s cos[2tt(VT) - <j> s ], (35) 
where R — lAU/c = 499.00478s and d0/dt is the azimuthal velocity of the orbit of the inspiralling object (see Eq. 

(DSD- 

IV. NOISE INDUCED PARAMETER ERRORS 

This section contains an overview of the basic elements of signal analysis that we will use in Section [V] to estimate 
the accuracies with which GW observations will be able to determine the system parameters. 

The measured strain s(t) in a GW detector is a time series that contains both a true GW signal, h(t), and 
instrumental noise, n(t). In the context of LISA, the output of the equivalent two arm Michelson detectors can be 
represented as 

s a (t) — h a (t) + n a (t), a = I, II. (36) 
Given two time series, we define the overlap as 



(p | q) = 2^ / \P*M)qM) +PMWM)] /Sn(f) df, (37) 



n 



where S n (f) stands for the one-sided spectral density of the instrumental noise. Assuming that each Fourier component 
of the noise, h a (f), is Gaussian distributed, and uncorrelated with other Fourier components (i.e., the noise is 
stationary) , the ensemble average of the Fourier components of the noise has the property 

(n a (f) Mf'T) = \&{f - f')SnU)5 a p. (38) 
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The former relation defines the spectral density S n (f), which is the same for the two Michelson data streams. Fur- 
thermore, the probability distribution function for the noise n(t) is given by 



/ (n |n )\ 

p(n = n ) oc exp I - I . (39) 

We can interpret Eq. (|39[1 as the probability that the actual noise realization is no. Assuming that we have made a 
detection, i.e., the output of the detector is given by s(t) = h(i; (9 truc ) + no(t), where no(t) is the specific realization of 
the noise and #true is the unknown true value of the parameters of the source, the likelihood of measuring the output 
signal is 



w olfl s ( (s-h(g truc )) I (s~h(e tIUC )) \ 

A (s 1 6»trueJ oc exp I I. (40) 

In order to reconstruct the most probable value of the parameters of the source, and compute their respective errors 
we need to calculate the probability of the parameters given the data, i.e., the posterior probability, which, according 
to Bayes's theorem, is given by the product of the likelihood function, Eq. (|40l) . and the prior probability po(^true), 
i.e., 



p (h(0 t rue) I s) = JVpo(0true)ejq> f (h(0 true ) I s) - - (h(0 true ) I h(S true ))l , (41) 

where the factor (s|s) /2) has been absorbed into the normalization factor N. For a given measured signal s(t), 
the gravitational waveform h(t) that best fits the data is the one that minimizes the quantity (s — h | s — h). This 
condition is also satisfied by the maximum likelihood estimators of the parameters, and corresponds to the point in 
parameter space with the highest SNR in a matched filtering search, namely, 



> fl,)l = TOT" < 42 » 

In the limit of high SNR, a locally flat prior would be a reasonable assumption because the best-fit parameters will 
have a Gaussian distribution centered on the correct values. Hence, we can expand Eq. (|41|) about the peak, 9 = #true, 
by setting 9 i = 8* + A9\ and find 



p(A9\s) = Me 



dh 



(43) 



where p(A9 \ s) is the Gaussian probability distribution of the parameter estimation errors A9, and is the Fisher 
Information Matrix. Additionally, for large SNR, the covariance of the posterior probability distribution, (r -1 ) 4 - 7 , 
gives the expectation value of the errors A9 1 



(A9 i A9 j ) = (T- 1 )^ + O(SNR)- 1 . (44) 

In our analysis, we will use a simplified definition of the inner product, Eq. (|37l) . For white noise, i.e., S n (f) = const., 
it takes the simple form 2S^ 1 J2 a I™ Pa(t)q a (t)dt, by Parseval's theorem [1§|- Following Barack & Cutler we 
can define the "noise-weighted" waveform as follows 

K{t) = ; f(t) = - — , (45 

7r at 



JS h (f(t)) 

and rewrite the Fisher matrix approximately as 



T 



I d a h a {t)d b h a (t)dt . (46) 
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A. Noise model 



The function Sh(f) is the total LISA noise, which has three components: instrumental noise, confusion noise from 
short-period galactic binaries, and confusion noise from extragalactic binaries. We use the same prescription as in 
Barack and Cutler [19], namely 

S h {f(t)) = min{5r t (/)/exp(- K T- i 1 ssion d7V/d/) + S^(J), 

Sr\f) + Sf(f) + S^ a \f)}, (47) 
where the various components have the following analytic forms [l9| 

S'h st (f) = 9.18 x l(T 52 /~ 4 + 1.59 x 1(T 41 + 9.18 x 1(T 38 / 2 Hz" 1 , 
Sf(f) = 2.1x10^ (X) Hz" 1 , 

£CX. gal = 4 2 x 1Q -47 ( _f_\ ' jj z -l (48) 

\ 1Hz / 

Here dN/df is the number density of galactic white dwarf binaries per unit GW frequency, and k is the average 
number of frequency bins that are lost when each galactic binary is fitted out. We use 



dN _ _ - 1 /1Hz 



df 

3 yr, and k — 4.5. 



11/3 



2xlO- 3 Hz" 1 ( — ) , K T-^ sion = 1.5/yr, (49) 



V. PARAMETER ESTIMATION ERROR RESULTS 



In this Section we explore the accuracy with which LISA observations will be able to measure the spin of the 
inspiralling object. To effectively address this problem, we will estimate noise-induced errors for a variety of binary 
systems using the formalism described in Sections |TT] and IIVI To estimate the noise-induced errors we will use the 
inverse Fisher Matrix (see Eq. (j44|)). We have verified that the Fisher Matrices used in our studies were convergent over 
several orders of magnitude in the offset used to compute the numerical waveform derivatives. The same convergent 
behaviour was exhibited by the associated inverse matrices. 

The noise-induced errors will be quoted for fixed values of the intrinsic parameters of the source, but with a Monte 
Carlo simulation over possible values of the extrinsic parameters. We compute the Fisher Matrix for a source at 
D = lGpc, and the corresponding SNR using the expression 




We then renormalise the results to a fixed SNR, which was chosen separately for each of the systems considered. This 
choice was based on a Monte Carlo simulation of the SNRs, as described in the next subsection. 



A. Determination of typical SNRs 



To estimate a suitable "typical" SNR for each type of source to use as a reference, we carried out a Monte 
Carlo simulation in which the extrinsic parameters of each source were chosen randomly, with the events distributed 
uniformly in comoving volume out to a redshift of z = 1, and detected in a certain time window at the detector, which 
we took to be one year. For each event, we computed the SNR and then looked at how the SNRs of the detected 
events were distributed. For systems involving the more massive inspiraling objects, events at redshifts z > 1 could 
still be detected by LISA. However, given the uncertainties that exist in the formation mechanisms for IMBHs and the 
inspiral rates, we chose to assume events that occur only for redshifts z < 1. This would be a conservative assumption 
if it were used to compute the event rate for detected inspirals. However, if there are a large number of inspirals 
occurring at higher redshift, the reference SNRs we compute in this way could be somewhat larger than those of the 
systems LISA eventually detects. A full assessment of the parameter estimation accuracies that LISA will achieve 
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would require detailed modelling of the source population that is beyond the scope of the current paper. The results 
presented here can be readily rescaled to other SNRs for any future study of this kind. 

In the Monte Carlo simulation, the intrinsic parameters of the source were fixed, with redshifted masses equal to 
the values quoted earlier, e.g., 1OM + 10 6 Af© etc. For each source, we computed the SNR accumulated in a one 
year observation prior to plunge. This was accomplished by choosing the initial radial coordinate, po, such that the 
inspiralling BH reached the last stable orbit after one year of inspiral. Figure [T] presents the normalized cumulative 
distribution functions of the SNR computed in this way, for binary systems with central BHs of redshifted mass 
lO 6 Af , and spin parameter q = 0.9. The inspiralling BHs have specific spin parameter x = 0.9, and redshifted 
masses fx = 1OM , fi = 1O 2 M , \x = 1O 3 M , and /i = 5x 10 3 M Q . In Fi gure[2j we present the normalized cumulative 
distribution function for binary systems with redshifted masses 5 x 1O 3 M + 10 6 M Q , for three additional combinations 
of the spin magnitudes of the binary components. 




SNR SNR 



FIG. 1. The normalized cumulative distribution function for the signal-to-noise ratio of a cosmological population of binary 
systems with central BHs of redshifted mass M = 1O 6 M0, and spin parameter q = 0.9. The inspiralling BHs have specific spin 
parameter x = 0.9, and redshifted masses /x = lOM© (top left panel), fi = 10 2 Mq (top right), fi — 1O 3 M0 (bottom left), and 
fi = 5 x 1O 3 M (bottom right). 



The way to interpret the normalised cumulative distributions of Figures [1] [5] is the following: the value of the 
cumulative distribution at a given SNR t , indicates the fraction of events with SNR < SNR t . Hence, for the least 
massive system, the cumulative distribution indicates that ~ 50% of events in the Monte Carlo sample have SNR <J 30. 
This "median" SNR is a suitable reference SNR at which to quote parameter estimation accuracies for the binary 
systems considered in our studies. 

Based on Figures [1] and [2] we will therefore normalise our results to the following reference SNRs: for the inspiralling 
BHs with masses fi — 1OM , 1O 2 M , lO 3 Af , 5 x 1O 3 M , and spin magnitude x — 0-9, we will use a reference SNR 
of 30, 150,400, 1000, respectively. For the more massive system, 5 x 1O 3 M + lO 6 Af , Figure [5] suggests that for the 
combination of spin magnitudes q = 0.9, x = 0.1, we should use a reference SNR of 1000; whereas for the combination 
q = 0.1, x = 0.9, and q = 0.1, x = 0.1, we should use SNR= 500. 

The results shown in Figures [1] and [5] provide information that could be used to renormalise the results presented 
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FIG. 2. Normalized cumulative distribution function for the signal-to-noise ratio of a cosmological population of binary 
systems with BHs of redshifted mass 5 x 10 3 M© + 10 6 Mq. The panels show results for three different combinations of the spin 
magnitudes of the central and inspiralling BHs, q, x, namely: top panel q = 0.9, x — 0.1; bottom-left panel q = 0.1, x = 0.9; 
bottom-right panel q — 0.1, x = 0.1. 

later in the paper, in Tables HHIVH to a cosmological population of sources at z <J 1. 

B. Parameter estimation results 



The parameter space we have considered is 11-dimensional. Five of these are intrinsic parameters, namely 
km, In M, q, x,Po> while the other six are extrinsic or phase parameters. We summarize the physical meaning 
of the parameters in Table fl] 

In Tables HH IIII1 IIV1 fVl and IVI1 we summarise the results of Monte Carlo simulations of the parameter estimation 
errors, over possible values for the extrinsic parameters of the source, but with fixed intrinsic parameters. Each Table 
present results for a different binary system. All of the systems have a central black hole with mass M — 10 6 M Q 
and spin parameter q — 0.9, and an inspiraling BH with spin parameter \ — 0.9, but we make four different choices 
for the mass of the inspiralling BH, /i = 10M Q , 10 2 M Q , 10 3 M Q , 5 x 10 3 M Q . For completeness in the analysis, and 
to explore the trend with the spin of the small/big body, we consider three additional cases for the most massive 
inspiraling BH, in which we vary the spin parameters of the components, choosing q — 0.9, \ = 0.1; q = 0.1, \ = 0.9; 
and g = x = 0.1. These results are shown in Tables W\ and IVll 

The results in Table |H] can be compared to out previous analysis for EMRIs that ignored the spin of the smaller 
black hole [lH. We see that i) including the spin of the small CO for EMRIs with mass ratios rj <J 10~ 5 will not 
significantly affect parameter determination or detection and, ii) GW observations will not be able to constrain at all 
the spin parameter of the inspiralling BH for systems with such a small mass ratio. These results are consistent with 
the arguments presented by Barack & Cutler [19| in the sense that including small body spin effects for EMRIs has 
a minor effect on the orbital evolution of the system. 
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In n 


mass of inspiralling object 


InM 


mass of central SMBH 


q 


magnitude of (specific) spin angular momentum of SMBH 


X 


magnitude of (specific) spin angular momentum of inspiralling object 


Po 


Initial radius of inspiralling object's orbit 


4>o 


Initial phase of inspiralling object's orbit 


6s 


source sky colatitude in an ecliptic-based system 


<t>s 


source sky azimuth in an ecliptic-based system 


9 K 


direction of SMBH spin (colatitude) 




direction of SMBH spin (azimuth) 


InD 


distance to source 



TABLE I. This table describes the meaning of the parameters used in our model. The angles (9 s, <j>s) and (9 k, 4>k) are defined in a fixed 
ecliptic— based coordinate system. 



Model 


Distribution of log 10 (AX) in error, AX, for parameter X — 


ln(m) 


ln(M) 


q 


X 


po 


4>o 


6s 


<t>s 


Ok 


4>K 


ln(D) 


q = 0.9 
X = 0.9 


Mean 


-3.99 


-3.58 


-4.13 


1.68 


-3.31 


-0.99 


-1.54 


-1.58 


-1.14 


-1.07 


-1.09 


St. Dev. 


0.110 


0.122 


0.114 


0.098 


0.124 


0.453 


0.133 


0.199 


0.433 


0.456 


0.311 


L. Qt. 


-4.08 


-3.61 


-4.34 


1.63 


-3.35 


-1.35 


-1.63 


-1.71 


-1.49 


-1.40 


-1.37 


Med. 


-3.99 


-3.55 


-4.08 


1.68 


-3.29 


-1.16 


-1.53 


-1.62 


-1.24 


-1.08 


-1.16 


U. Qt. 


-3.89 


-3.51 


-3.95 


1.74 


-3.23 


-0.85 


-1.49 


-1.51 


-0.85 


-0.79 


-0.88 



TABLE II. Summary of results of the Monte Carlo simulation of Fisher Matrix errors for spinning BH systems with specific 
spin parameters q = x = 0.9, and masses fj, = IOMq, M = 10 6 A/q. We show the mean, standard deviation, median and 
quartiles of the distribution of the logarithm to base ten of the error in each parameter. Results are given for the kludge model 
with conservative corrections to 2PN order. The angles <f>o and ao, specifying LISA's position and orientation at t = 0, are set 
to zero. Note that the results have been normalised to fixed SNR= 30. 

For systems with mass ratios rj ~ 10 -4 , Table [TTTI indicates that it will not be possible to measure the spin of 
inspiralling BHs in this regime either. However, Table IIVI indicates that it should be possible to measure the spin 
of the inspiralling BH for systems with mass ratios r\ ~ 10~ 3 . At an SNR of 400, GW observations should be able 
to determine the spin parameter of the inspiralling BH for systems with component masses 10 3 M Q + 10 6 M Q to an 
accuracy of ~ 30%. 

For inspiralling BHs of mass /i = 5x 10 3 Af Q , the results are quite similar, but the higher expected SNR for such 
systems should allow us to measure the small object spin to higher precision. For these systems, Table IVl and Figure[3] 
indicate that, at a fixed SNR= 1000, LISA measurements will be able to determine the inspiraling BH mass, SMBH 
mass, SMBH spin parameter, and inspiraling BH spin parameter with fractional errors of ~ 10~ 3 , 10~ 3 , 10~ 4 , and ~ 
10 _1 . We also expect to determine the location of the source in the sky and the SMBH spin orientation to ~ 10~ 4 
steradians. 

We conclude that the inclusion of small body spin effects could become relevant for binaries with mass ratios 
rj ;> 10~ 3 . Having found this threshold, we can explore the trend with the spin of the small/big body in this case. 
Tables IVl and PVTl present four different combinations of the spin parameters of the more massive binary's components. 
For q = x = 0.9, Table IVl suggests that GW observations will be able to measure the spin of the inspiraling body to a 



Model 


Distribution of log 10 (AA) in error, AX, for parameter X — 


ln(m) 


ln(M) 


1 


X 


Po 


4>o 


e s 




6 K 


4>K 


HP) 


q = 0.9 
X = 0.9 


Mean 


-3.78 


-3.62 


-4.98 


0.19 


-2.76 


-1.60 


-2.04 


-2.00 


-1.81 


-1.68 


-1.78 


St. Dev. 


0.074 


0.075 


0.112 


0.068 


0.075 


0.330 


0.206 


0.246 


0.361 


0.396 


0.270 


L. Qt. 


-3.84 


-3.68 


-5.06 


0.14 


-2.82 


-1.86 


-2.20 


-2.18 


-2.08 


-1.98 


-2.02 


Med. 


-3.79 


-3.63 


-4.98 


0.19 


-2.77 


-1.66 


-2.03 


-2.01 


-1.89 


-1.74 


-1.85 


U. Qt. 


-3.74 


-3.57 


-4.88 


0.25 


-2.71 


-1.43 


-1.87 


-1.88 


-1.56 


-1.40 


-1.58 



TABLE III. As Table UTJ but for an inspiralling BH with mass /U = 100M Q . Results are quoted at a fixed SNR of 150. 
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Distribution of log 10 (AX) in error, AX, for parameter X = 



Model 


ln(m) 


ln(M) 


q 


X 


Po 


4>o 


e s 


0S 


9 K 


4>K 


ln(D) 


q = 0.9 
X = 0.9 


Mean 


-3.33 


-3.15 


-4.62 


-0.52 


-2.03 


-1.39 


-1.88 


-1.84 


-1.67 


-1.56 


-1.73 


St. Dev. 


0.096 


0.097 


0.091 


0.084 


0.097 


0.407 


0.344 


0.363 


0.444 


0.470 


0.321 


L. Qt. 


-3.41 


-3.23 


-4.70 


-0.59 


-2.11 


-1.66 


-2.15 


-2.06 


-1.95 


-1.87 


-1.94 


Med. 


-3.38 


-3.20 


-4.67 


-0.55 


-2.08 


-1.50 


-1.82 


-1.79 


-1.72 


-1.62 


-1.82 


U. Qt. 


-3.24 


-3.06 


-4.54 


-0.45 


-1.93 


-1.26 


-1.60 


-1.68 


-1.43 


-1.30 


-1.55 



TABLE IV. As Table HU but for an inspiralling BH with mass fi = 10 3 M Q . Results are quoted at a fixed SNR of 400. 



Model 


Distribution of log 10 (AX) in error, AX, for parameter X = 


ln(m) 


ln(M) 


1 


X 


Po 


4>o 


e s 


0S 


6k 




HP) 


q = 0.9 
X = 0.9 


Mean 


-3.12 


-2.94 


-4.38 


-1.07 


-1.60 


-1.48 


-2.06 


-2.07 


-1.85 


-1.71 


-1.90 


St. Dev. 


0.089 


0.089 


0.083 


0.081 


0.089 


0.417 


0.379 


0.391 


0.434 


0.468 


0.326 


L. Qt. 


-3.14 


-2.99 


-4.40 


-1.09 


-1.65 


-1.91 


-2.39 


-2.23 


-2.13 


-2.00 


-2.10 


Med. 


-3.13 


-2.94 


-4.38 


-1.06 


-1.64 


-1.72 


-2.07 


-2.02 


-1.89 


-1.73 


-1.95 


U. Qt. 


-3.03 


-2.84 


-4.37 


-1.02 


-1.54 


-1.42 


-1.87 


-1.74 


-1.62 


-1.53 


-1.71 


q = 0.9 

x = o.i 


Mean 


-2.62 


-2.44 


-3.92 


-0.60 


-1.19 


-1.49 


-2.09 


-2.04 


-1.75 


-1.62 


-1.81 


St. Dev. 


0.181 


0.185 


0.176 


0.172 


0.205 


0.359 


0.440 


0.365 


0.398 


0.444 


0.411 


L. Qt. 


-2.77 


-2.59 


-4.07 


-0.70 


-1.31 


-1.96 


-2.44 


-2.32 


-2.12 


-2.09 


-2.10 


Med. 


-2.60 


-2.41 


-3.90 


-0.61 


-1.21 


-1.66 


-2.06 


-2.03 


-1.82 


-1.77 


-1.89 


U. Qt. 


-2.46 


-2.29 


-3.77 


-0.44 


-1.09 


-1.36 


-1.78 


-1.69 


-1.60 


-1.41 


-1.67 



TABLE V. As Table HU but for an inspiralling BH with mass fx = 5 X 10 3 M Q . Results are quoted at a fixed SNR of 1000. 

precision slightly better than ten percent. For q = 0.9, \ = 0.1, this same Table shows that when the inspiraling body 
is slowly rotating, the accuracy with which we can determine the intrinsic parameters of the system is a factor of ~ 3 
times worse. However, even in this latter case, the small body spin magnitude could be measured to a fractional error 
of 25% and the accuracy with which the other extrinsic parameters can be determined remains basically unchanged. 

Table PVTl shows that when the central SMBH is slowly rotating i) the accuracy with which the intrinsic parameters 
of the system can be determined is not very sensitive to the spin of the inspiraling body; ii) GW observations will not 
provide an accurate measurement of the spin of the inspiraling object. These results are to be expected — the influence 
of spin couplings on the orbital evolution is what allows the spin of the inspiraling object to be measured from the 
gravitational wave emission. These spin couplings are enhanced when the two objects are more rapidly rotating, and 
when the inspiraling object is in the strong-field region close to the central massive object. As the innermost stable 
orbit is closer to the central SMBH when the latter is more rapidly rotating, this further increases the measurability 
of the spin coupling in that case. We conclude that the inspiraling object's spin will be most measurable when both 



Model 


Distribution of log 10 (AX) in error, AX, for parameter X = 


ln(m) 


ln(M) 


1 


X 


po 


00 


6s 


0s 


9k 




HP) 


q = 0.1 
X = 0.9 


Mean 


-3.09 


-2.92 


-2.61 


-0.15 


-1.62 


-1.58 


-1.99 


-1.89 


-1.75 


-1.63 


-1.79 


St. Dev. 


0.071 


0.072 


0.064 


0.062 


0.072 


0.365 


0.401 


0.416 


0.387 


0.416 


0.283 


L. Qt. 


-3.16 


-2.99 


-2.66 


-0.20 


-1.68 


-1.82 


-2.25 


-2.15 


-2.00 


-1.92 


-1.98 


Med. 


-3.09 


-2.93 


-2.62 


-0.16 


-1.62 


-1.66 


-1.87 


-1.85 


-1.77 


-1.64 


-1.83 


U. Qt. 


-3.01 


-2.85 


-2.56 


-0.11 


-1.56 


-1.41 


-1.67 


-1.64 


-1.54 


-1.38 


-1.62 


q = 0.1 
X = 0.1 


Mean 


-3.09 


-2.92 


-2.61 


-0.15 


-1.62 


-1.68 


-1.99 


-1.90 


-1.85 


-1.71 


-1.86 


St. Dev. 


0.069 


0.070 


0.064 


0.060 


0.070 


0.282 


0.397 


0.444 


0.348 
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U. Qt. 
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-1.61 
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-1.72 



TABLE VI. As Table fVl but for a slowly rotating central BH with spin parameter q = 0.1. Results are quoted at a fixed SNR 
of 500. 
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FIG. 3. Distribution of parameter measurement error estimates, A, computed in the Monte Carlo simulations. We show results 
for the system with (i = 5x 10 3 Mq, M = 1O 6 M0, q — 0.9, \ ~ 0-9- The panels show, from left to right, the error distributions 
for, top row: A(log 10 (ln/i)), A(log 10 (ln M)); and bottom row: A(log 10 (g)) and A(log 10 (x)) respectively. Results are quoted at 
fixed SNR of 1000. 



bodies are rapidly rotating, which is supported by the results in Tables [V] and EH We see that GW observations 
in the low-frequency band could provide useful information not only on the spin distribution of SMBHs, but also of 
intermediate mass BHs that exist in the centres of galaxies, if binaries comprising two rapidly spinning objects are 
observed with sufficient SNR. 

One final thing to note from Tables Hill VII is that the precision with which the other intrinsic parameters can be 
measured decreases as the mass of the inspiraling object increases. This is true even though the reference SNR used for 
each table is increasing, and would be even more pronounced if the results were renormalised to the same reference SNR 
for all systems. This trend is due to the fact that the more massive the binary system is, the more rapid the inspiral 
and the less time the inspiralling object spends close to the central SMBH. For the system fj, = WMq + 10 6 M Q with 
q = x = 0.9, the smaller object spends the last six months of inspiral within the region p € (~ 6.01A/ — >~ 2.32M). 
In contrast, for the system fj, = 5 X 1O 3 M + 10 6 M Q with q = \ — 0-9) the object is at p s» 29.5M six months 
before plunge and hence spends much less time close to the innermost stable circular orbit (ISCO) of the central BH. 
Furthermore, the number of cycles completed over the last year of inspiral by the latter system is only ~ 10% of the 
number of cycles completed by the former. 
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VI. MODEL-INDUCED PARAMETER ERRORS 



The Fisher Matrix formalism we described in Section IIVI allowed us to estimate parameter estimation errors that 
would arise from noise in the detector. We shall now explore errors that would arise from using an approximate 
waveform model for data analysis. Even in the absence of noise, for a particular true signal, s(t), the kludge waveform 
model, h(t), that is the best match may have different parameters to the true waveform, which introduces another 
parameter error. These are the "model" errors. Using the framework developed by Cutler and Vallisneri [2l|, and 



our own analysis on model errors for EMRIs in the context of spinless particles 18(, we will present results on the 
magnitude of the model errors that could arise in E/IMRIs in which the inspiral component has significant spin. 

As discussed by Cutler & Vallisneri [2l|, we have no way to know what the "true" signal is. However, we can use 
their formalism to estimate how important for parameter estimation the inclusion of various terms in the waveform 
model may be. In this section, we look specifically at what effect ignoring various "conservative" post-Newtonian 
terms will have on parameter determination. One feature that distinguishes "model" errors from noise-induced errors 
is that whereas the former are SNR independent, the latter scale as 1/SNR (see (|44)l ). Hence, it might be the case that 
for the systems with highest SNR, theoretical errors could dominate the total parameter-estimation error. In order 
to find out whether including conservative corrections is important for detection and parameter estimation, we will 
build waveform templates that include all, only part, or none of the conservative pieces we derived in Section [TTJ The 
kludge waveform template that includes all terms at 2PN order, see Eq. (|28p. will be taken as the "true" waveforms 
/igr. This model includes all the coefficients do,di..., etc., of the azimuthal frequency given in Eq. ([28"]). We then 
estimate the model errors by searching for /igr using templates, /iap, that include none or only part of these terms. 
When we refer to a template including conservative corrections at 1.5PN order, we will mean a template in which 
the angular frequency is given by Eq. (|28l) including the coefficients do, d\, di.5, /1.5 (note that we found 51.5 = 0). 
Similarly a template including conservative corrections at OPN order, will mean a template in which all the coefficients 
except the leading 1 inside the curly bracket of Eq. (|2"8l are set to zero. We refer to this as testing the importance of 
including conservative corrections, as we will only alter the terms entering expression (|28l) for the frequency, and will 
not alter the radiative part of the waveform model described by the radial evolution in Eq. (IT2l) - (fT5|) . The reason for 
treating these two parts of the model differently is that the latter, radiative, part of the evolution can be computed 
straightforwardly from solution of the Teukolsky equation, while a calculation of the conservative part of the evolution 
requires knowledge of the local self-force, which is much more difficult to compute and is presently unknown. 

At the end of this section, we will also assess the overall importance of including the spin of the small body in the 
model by comparing a waveform with x — 0, to one in which \ is kept as a free parameter. This will directly test 
whether using templates that do not include the small body spin will be adequate for accurate parameter estimation 
from GW observations. The rule of thumb in both of these analyses will be that if omitting part of the model gives 
rise to a parameter error that is comparable to or smaller than the noise-induced error, then we can safely ignore that 
part of the model when constructing search templates. 

In order to provide an analysis as self-contained as possible, we shall now present a brief overview of the Cutler & 
Vallisneri 21^ model error formalism that will be used for our subsequent analysis. Consider two manifolds embedded 
in the vector spaces of data streams. One of them is covered by the true waveforms {/igr(# 1 )} and the other one by 
the approximate waveforms {/iap(0 z )}- Given a signal s = hcR^) + n, the best fit 9 l is determined by the condition 



(^•h AP (^ 



s-h AP (0 l )J=O. (51) 
Furthermore, at first order in the error A0 l (n) = 8 l (n) — 9\ Eq. (IBTj) takes the form 



r- 1 ^))^ (^h AP (0)| n) + (r- 1 ^)) 1 ' (0,-h A p(0)| h GR (0) - h AP (0)) , (52) 



where the Fisher matrix is evaluated using the approximate waveforms Tij(6) = (9ih A p(0)|c?jh A p(0)). 

Relation (f52"|) indicates that, at leading order, A9 l is the sum of two contributions. The first one is due to noise in 
the detector, A n 6 l , whereas the second one, A t h# 4 , is the contribution due to the inaccurate waveform model. These 
are given, respectively, by 

A„0* = (r-^y'^hAP^In), A th 9 l = (r- 1 (0)) iJ (5 J h AP (0)|h GR (0)-h A p(0)). (53) 



If we knew both 9 and the noise realization n, then these equations would allow us to determine 9. However, 
experimentally we are only able to determine the h A p(#) that is the best fit waveform for the data stream, s, and we 
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are unsure about the error A6 = 6 — 9. In addition, we do not know 9 in Eq. (|53|) . At leading order, we can replace, 
^Gk{9) - h A p(0) by h GR (0) - h A p(6>), obtaining 



A th l = (r- 1 (0)) U (a> AP (0)|h GR (0) - h AP (9) 



(54) 



This relation is both noise and SNR independent. This property, along with the fact that A t h$ 4 is not averaged out if 
the same event is measured by a large number of nearly identical detectors leads us to consider Ath$ l as a systematic 
error. 

Cutler and Vallisneri [2l[ found that the leading order approximation, Eq. (j54|) . was not very good, unless the 
waveform is re-written in an amplitude-phase form 



h a {f) = A Q (/)e 4 *" (/) 



The amplitude A and phase "J are given by 



Ai = ^Ja 2 + F 2 + +A 2 x F 2 x , *z = 20 + Vz, 



2D 



A TJ = -^LJMF; 



2D 



ii. 



A 2 F 2 

^■x r II,x ' 



*// = 2<j> + ipu, 



where A +)X are given by (j3Tj) with A + = A x = 2p 2 <ft 2 , i ? {/j/ ; + ! x} are given by (|32|) and (|33l) . and 



ipi 



arctan 



A x Fr 



A+Fj. + 

The first order approximation to this expression 



ipil = arctan 



A x Fn,x 
A + F n ,+ 



Vtht 



r- 1 ^))'^ [AA + iAA*]e l * a,hAP(0) 



(55) 

(56) 
(57) 

(58) 
(59) 



at 9 



was found to give reliable results when compared to more accurately computed error estimates [2l| , so we use this form 
again here. Equation (|59"j) behaves better than Eq. (f54"|) since the difference between two waveforms, 1iap(0) — 1iap(#), 
is not very well approximated by the first term in its Taylor expansion. The differences in both the amplitude and 
phase of the waveform are individually well approximated by the linear terms in the Taylor series [21] . In fact, 
(1541 is reliable only as long as the phase difference between the two waveforms is much less than one radian, i.e. 
A^ <9j \&ap (#) <C 1, whereas for (|59"|) we just require A9 l A9^ didj^ ap(9) <C 1. This condition is much less restrictive 
than the former one. 

We can now use equation (I59|) to estimate the magnitude of the parameter errors that arise from inaccuracies in 
the template waveform. In j 1 81 ] we introduced a waveform model that included conservative self-force corrections at 
2PN for Kerr circular-equatorial EMRIs. The conservative corrections used in that model were derived using the 
same method employed in Section HH since accurate fully relativistic corrections are not currently available for objects 
moving in a Kerr background. There is, however, an ongoing program that aims at developing accurate and efficient 
numerical methods for computing such corrections. As mentioned in Section [Vj first-order self-force results have 
been derived for a scalar charge moving on circular and eccentric equatorial orbits in the Kerr spacetime [III [2i| . 
The extension of this analysis to the gravitational case is now under investigation. The computation of conservative 
corrections for the case we consider in this paper, namely, spinning BHs moving in the background of a Kerr BH, is 
beyond the scope of the current self-force program, which is focused on inspiraling objects that are non-spinning. 

In the absence of fully perturbative calculations, we can use the kludge model to assess the importance of conservative 
corrections for detection and parameter estimation. This investigation will indicate whether it is necessary, for source 
detection and parameter estimation, to know these corrections accurately before we can properly analyse data from 
a space-based gravitational wave detector. This is an important question considering that these corrections are not 
yet known. 

We will compute the model error that arises when we omit some or all of the "conservative" corrections in the 
waveform template and hence calculate the ratio 1Z of this model error to the error that will arise from noise in the 
detector. This ratio will indicate the importance of including the conservative corrections for parameter determination. 
If 1Z <; 1, then parameters estimated using a model that ignores the conservative corrections should still be reliable, 
but if 1Z » 1 then it is clear that we must include the conservative corrections. The ratios of parameter errors 
to Fisher Matrix errors obtained from the Monte Carlo simulations over extrinsic parameters are summarized in 
Tables EED rvTnl El H EH and EH 
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We quote results for the same four test systems that we used previously, and consider two different comparisons, 
namely, we take the "true" waveform to be our kludge waveform with 2PN conservative corrections, and the template 
to be a kludge waveform with either no conservative corrections ("OPN") or with conservative corrections to 1.5PN 
order ("1.5PN"), and do pairwise comparisons. The exact meaning of OPN etc. was described in detail at the start 
of this section. 
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TABLE VII. Summary of Monte Carlo simulation results for the ratio of model errors to noise-induced errors, computed using 
the Fisher Matrix, for spinning BH systems with n — IOMq. We show the mean, standard deviation, median and quartiles 
of the distribution of the logarithm to base ten of the ratio for each parameter. Results are given for various comparisons, as 
indicated and described in the text. A comparison "A vs B" uses model A as the true waveform, and model B as the search 
template. Note that the noise-induced errors are quoted at a fixed SNR= 30. 
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TABLE VIII. As Table IVII1 but for a spinning BH with mass fi = lOOAf© . Noise-induced errors are quoted at a fixed 
SNR= 150. 

Tables IVIIIIXII indicate that, as expected, the ratio 1Z of model errors to noise-induced errors becomes smaller as 
the approximate waveform, Iiap, becomes closer to the "true" waveform, hop, . The improvement going from OPN to 
1.5PN is largest for the more massive binaries. 

From Table IVTll we see that for the binary systems lOM© + 10 6 Mq, the vast majority of sources fulfill the condition 
that the model errors are smaller than the noise-induced errors, for the particular SNR, 30, at which these results are 
quoted. We note also that the error ratio for the spin of the inspiralling object appears to be similar to the ratio for 
the other intrinsic parameters. However, this value is the ratio of two large numbers. As indicated in Section [Vl the 
median of its noise-induced error is of the order of ~ 10 1 ' 6 , whereas its model error is of the order of ~ 10 1,2 . Hence, 
although the ratio looks small, this is not an indication that this parameter can be accurately measured by EMRI 
observations. 

Tables l\THT[Xll indicate that, for the more massive BH binaries, model errors are likely to be larger than the 
statistical errors. This result was already pointed out by Cutler & Vallisneri [2l| in the context of the inspiral of 
non-spinning massive BH binaries. In their studies, they used a simple PN model with no spin corrections, for 
sources with a duration of 1 year, which were truncated at the ISCO of non-spinning BHs, i.e., p = 6M. We have 
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TABLE IX. As Table lVlfl but for a spinning BH with mass /i = 1O 3 M . Noise-induced errors are quoted at a fixed SNR= 400. 
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TABLE X. As Table IVTT1 but for a spinning BH with mass // = 5 X 1O 3 M0. Noise-induced errors are quoted at a fixed 
SNR— 1000. 



now extended that initial calculation by i) building a more accurate waveform model using the equations of motion 
derived by Saijo et al. [l| for spinning inspiralling objects; ii) augmenting this model with spin-orbit and spin-spin 
couplings taken from perturbative and PN calculations to amend the equations of motion; and iii) modelling the 
inspiral phase using the most accurate fluxes available which include perturbative corrections for the spin of the 
inspiralling body, spin-spin couplings and higher-order fits to Teukolsky calculations. Additionally, we have used a 
consistent prescription for the computation of the ISCO at which we truncate the waveforms. Small body spin effects 
produce fairly negligible changes in the value of ISCO for IMRIs, but we have included these in order to be consistent 
throughout our analysis. 

There are two reasons why the model error to noise-error ratio appears to increase as the mass of the inspiralling 
body increases. Firstly, we are using a higher reference SNR to quote the noise errors for the more massive systems, 
since these are intrinsically louder. At higher SNR, the noise-induced error decreases while the model error is fixed. 
Nonetheless, if we renormalise the results in the tables to a fixed SNR, across all sources, we still see that the error 
ratio, TZ, for the 2PN to OPN comparison, tends to increase as the mass ratio rj is increased. This is to be expected, 
since the corrections that we are omitting in the comparison are proportional to the mass-ratio and therefore should 
have a greater impact for higher mass-ratio systems. 

However, Tables IVIIIIiXll also show that for the 2PN vs 1.5PN comparisons, the error ratio, TZ, increases, but 
remains small and if the results were renormalised to a fixed SNR, it would actually decrease slightly as the mass of 
the inspiraling object was increased. For the 2PN vs 1.5PN comparison with \ = 0.9 and // = 5 x 1O 3 M0, the error 
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Model 


logic °f the ratio 1Z of model to noise-induced error for parameter X = 


ln(m) 


ln(M) 


<7 


X 


Po 


00 


0s 


0S 


6 K 




ln(D) 


2PN vs 1.5PN 
X = 0.9 


Mean 


0.41 


0.40 


0.49 


0.49 


0.38 


0.79 


0.66 


0.62 


0.88 


0.91 


0.87 


St. Dev. 


0.597 


0.615 


0.517 


0.643 


0.607 


0.698 


0.739 


0.723 


0.776 


0.801 


0.799 


L. Qt. 


0.05 


0.05 


0.15 


0.13 


0.02 


0.22 


0.17 


0.09 


0.16 


0.21 


0.12 


Med. 


0.46 


0.46 


0.60 


0.58 


0.45 


0.81 


0.73 


0.69 


0.93 


0.92 


0.99 


U. Qt. 


0.78 


0.75 


0.90 


0.87 


0.75 


1.47 


1.23 


1.19 


1.81 


1.75 


1.74 


2PN vs 1.5PN 
X = 0.1 


Mean 


0.43 


0.41 


0.51 


0.50 


0.41 


0.78 


0.66 


0.64 


0.90 


0.93 


0.92 


St. Dev. 


0.624 


0.619 


0.632 


0.681 


0.609 


0.702 


0.757 


0.695 


0.758 


0.749 


0.784 


L. Qt. 


0.06 


0.03 


0.17 


0.17 


0.03 


0.18 


0.20 


0.15 


0.18 


0.22 


0.19 


Med. 


0.51 


0.48 


0.63 


0.64 


0.46 


0.79 


0.74 


0.74 


0.94 


1.00 


0.96 


U. Qt. 


0.81 


0.83 


0.95 


0.93 


0.83 


1.49 


1.28 


1.24 


1.79 


1.83 


1.85 



TABLE XI. As Table E but for a central SMBH with spin parameter q — 0.1. Noise-induced errors are quoted at a fixed 
SNR= 500. 



log 10 of the ratio 1Z of model to noise-induced error for parameter X — 



Model 


ln(m) 


ln(M) 


q 


X 


Po 


00 


0s 


4>s 


0K 


4>K 


ln(D) 




Mean 


0.64 


1.12 


0.87 


6.21 


1.88 


3.41 


3.20 


3.21 


3.72 


3.98 


3.79 


X = 0.9 


St. Dev. 


0.519 


0.569 


0.610 


0.603 


0.629 


0.725 


0.617 


0.574 


0.697 


0.621 


0.649 


vs 


L. Qt. 


0.37 


0.99 


0.63 


5.94 


1.70 


3.33 


3.00 


2.99 


3.31 


3.50 


3.34 


x = o 


Med. 


0.73 


1.25 


1.02 


6.29 


1.98 


3.35 


3.29 


3.25 


3.81 


3.99 


3.85 




U. Qt. 


0.98 


1.43 


1.29 


6.58 


2.21 


3.55 


3.41 


3.44 


4.16 


4.33 


4.39 



TABLE XII. Model errors that arise from omitting the spin of the inspiralling object for systems with component masses 
[5 x 10 3 + 1O 6 ]M . 



ratio 7Z <J 4 in all cases at a reference SNR of 1000, which is already a manageably small value. At an SNR of 250, 
the error ratio would be less than 1. 

In this same case, the value of 1Z is a factor of ~ 20/30 larger for the intrinsic/extrinsic parameters for the 2PN 
vs 0PN comparison relative to the 2PN vs 1.5PN comparison. Altogether, these results suggest that a model which 
included conservative corrections up to 2PN order might be sufficient for parameter estimation, since the relative 
importance of the 1.5PN — > 2PN change is small, even for the most massive systems. More work is required to 
confirm this, by comparing the 2PN model to a higher order (3PN or 3.5PN) model, but the results presented here 
seem promising. 

In Table IXIII we show the model errors that would arise if we did not include the spin of the inspiralling body in 
the waveform template, but it was included in the "true" waveform. We show results for the system with component 
masses [5 x 10 3 + 10 6 ]M Q since the spin has the biggest impact in that case. This exercise is important in order to 
understand how much the parameter estimation we could achieve from a GW observation would be degraded if we 
did not use spinning templates in our search. This table shows that not including small body spin corrections could 
significantly degrade the accuracy with which the CO mass, SMBH mass and SMBH spin may be determined, since 
the model errors associated with these parameters are a factor of ~ 5, 20, 10, bigger than the noise-induced errors, 
respectively. We note that the large quoted error associated with the small body spin parameter Xj is merely an 
indication that we cannot determine that parameter using a waveform template that does not include it in the first 
place. These results suggest that the small body spin corrections in IMRIs should not be ignored in data analysis, 
since doing so may significantly degrade our determination of the system parameters. 

In summary, this analysis shows that model errors could be a limiting factor in determining the system parameters, 
but only for IMRIs. Similarly, the small object spin should not be ignored in IMRI waveforms. However, constructing 
a waveform template that included conservative corrections up to 2PN order could be sufficient to reduce systematic 
errors to an acceptable level. Such templates should be able to constrain a source to a sufficiently small region of the 
parameter space that a follow-up using more accurate and computationally expensive waveform templates would be 
possible, if such templates are available. 
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VII. CONCLUSIONS 



In this paper we have developed a kludge waveform model that includes small body spin corrections in order to 
explore the ability of a future low-frequency gravitational wave detector, such as LISA, to measure the spin of BHs 
inspiralling into much more massive Kerr BHs. This model uses the equations of motion derived by Saijo et al. 

and includes first-order conservative self-force corrections to compute the evolution of the inspiralling object's 
orbital frequency. The trajectory of the inspiralling object is computed using prescriptions for the fluxes of energy 
and angular momentum [l7[, augmented with perturbative results including spin-orbit and spin-spin couplings [23] . 

This analysis has demonstrated that LISA observations will not be able to measure the spin of stellar mass COs 
inspiralling into SMBHs. This result is in accord with arguments presented elsewhere [l9j . However, including small 
body spin effects could be relevant for signal detection and parameter estimation for binaries with mass ratio 77 ^> 10 . 
At a fixed SNR of 1000, a LISA observation of a binary with masses 5 x 10 3 A/ Q + 10 6 Af Q whose components have 
specific spin parameter q — \ — 0-9, will be able to determine the CO and SMBH masses, the SMBH spin magnitude 
and the inspiralling BH spin magnitude, %, to within fractional errors of ~ 10~ 3 , 10~ 3 , 10~ 4 , ~ 10%, respectively. 
We also expect to determine the location of the source in the sky and the SMBH spin orientation to within ~ 10~ 4 
steradians. Small body spin effects will be measurable when both components of the binary are rapidly rotating, but 
these effects will not be measurable or have a significant impact on parameter estimation when the central SMBH is 
slowly rotating. 

We have also studied in detail the importance, for parameter estimation, of including corrections that arise from 



the "conservative" part of the self-force, using the formalism developed by Cutler & Vallisneri [21] . We have found 
that, for a source with component masses 5 x 10 3 M Q + 1O 6 M0, and spin magnitudes q — x = 0-9, the relative 
error when using a 1.5PN template to detect a 2PN signal is small at fixed SNR, for all the systems we considered. 
Indeed, for the 2PN vs 1.5PN comparison, the ratio 7Z of the model errors to the noise-induced errors is of order 
1Z <; 4 for the intrinsic parameters, even for the "worst case" of a 5 x 1O 3 M + 1O 6 M binary detected with SNR 
of 1000. In contrast, when comparing 0PN to 2PN templates, this error ratio 1Z increases by a factor of ~ 20/30 
for intrinsic/extrinsic parameters. This suggests that including these corrections up to 2PN order might already be 
sufficient to reduce these systematic errors to an acceptable level. We also investigated the importance of including 
the spin of the small body in the waveform model. For systems with the most massive inspiraling objects, the model 
error that arises from omitting the spin could be one to two orders of magnitude larger than the noise-induced error, 
suggesting that including the spin in waveform templates will be important for such systems. 



These results extend an analysis initially carried out by Cutler & Vallisneri 21[ . Our studies include various 
ingredients that were not considered in there analysis. They used a purely post-Newtonian waveform, and compared 
a consistent 3PN approximate waveform to a "true" 3.5PN waveform. Moreover, their analysis included spin-orbit 
coupling as a model parameter, but the "true" waveforms were taken to be non-spinning in all cases, and they mostly 
focussed on near equal-mass systems. In our work, by contrast, we have focussed on more extreme-mass-ratio systems, 
we have included spin for the "true" waveforms, and our waveform model is based on accurate equations of motion 
for circular-equatorial spinning BH binary systems, which means it can be more readily applied in the strong-field. 
Finally, we have only considered the omission of "conservative" corrections when building our approximate templates, 
as these these are much more difficult to compute than the radiative corrections in the extreme-mass-ratio regime. 

The 5 x 10 3 Af Q + 1O 6 M system described here is an intermediate-mass-ratio inspiral. Another type of IMRI, in 
which a stellar-mass object (~ 1-M©) falls into an IMBH (~ IOOOMq), could be detected by future ground-based 
detectors, such as the Einstein Telescope (33|. In [34], HH, we assessed the capability of the Einstein Telescope 
to determine the parameters of such systems, when the inspiralling BH was non-spinning. The small-body-spin 
corrections in the model described here could be used to augment the IMRI model used in [34| to assess the importance 
and measurability of the small body spin effects in such systems. This should be studied in the future. 
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